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Abstract 


Pre-exposure prophylaxis (PrEP) is an important pillar to prevent HIV transmission. 
Because of experimental and clinical shortcomings, mathematical models that integrate 
pharmacological, viral- and host factors are frequently used to quantify clinical efficacy of 
PrEP. Stochastic simulations of these models provides sample statistics from which the clin- 
ical efficacy is approximated. However, many stochastic simulations are needed to reduce 
the associated sampling error. To remedy the shortcomings of stochastic simulation, we 
developed a numerical method that allows predicting the efficacy of arbitrary prophylactic 
regimen directly from a viral dynamics model, without sampling. We apply the method to var- 
ious hypothetical dolutegravir (DTG) prophylaxis scenarios. The approach is verified against 
state-of-the-art stochastic simulation. While the method is more accurate than stochastic 
simulation, it is superior in terms of computational performance. For example, a continuous 
6-month prophylactic profile is computed within a few seconds on a laptop computer. The 
method’s computational performance, therefore, substantially expands the horizon of feasi- 
ble analysis in the context of PrEP, and possibly other applications. 


Author summary 


Pre-exposure prophylaxis (PrEP) is an important tool to prevent HIV transmission. How- 
ever, experimental identification of parameters that determine prophylactic efficacy is 
extremely difficult. Clues about these parameters could prove essential for the design of 
next-generation PrEP compounds. Integrative mathematical models can fill this void: 
Based on stochastic simulation, a sample statistic can be generated, from which the pro- 
phylactic efficacy is estimated. However, for this sample statistic to be accurate, many sim- 
ulations need to be performed. 

Here, we introduce a numerical method to directly compute the prophylactic efficacy 
from a viral dynamics model, without the need for sampling. Based on several examples 
with dolutegravir (DTG) -based short- and long-term PrEP, as well as post-exposure pro- 
phylaxis we demonstrate the correctness of the new method and its outstanding computa- 
tional performance. Due to the method’s computational performance, a number of 
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Introduction 


Since its transfer to human in the early 20th century [1], HIV remains a major public health 
threat. According to UNAIDS estimates, approximately 38 million individuals worldwide are 
infected with the human immunodeficiency virus (HIV) [2]. HIV continues to spread and the 
latest incidence estimates amount to about 1.7 million new infections in 2019 [2]. Sub-Saharan 
Africa is hit hardest by the HIV pandemic, and due to COVID many services, including HIV 
control and treatment, had been suspended, which could lead to a long-term re-surge in infec- 
tions [3]. 

Nowadays, about 30 antiviral compounds are available that can stop HIV replication and 
prevent the acquired immunodeficiency symdrome (AIDS) and AIDS-related death [4]. How- 
ever, unlike many other infections, no cure is available to clear HIV, which can persist in latent 
reservoirs for decades [5, 6]. Available treatments therefore have to be taken life-long to pre- 
vent the relapse of virus from latent reservoirs and to prevent AIDS. As a consequence, much 
focus around fighting HIV turned towards HIV prevention. While a vaccine would be the 
ideal tool for the purpose, intrinsic difficulties have so far precluded the development of an 
effective vaccine against HIV [7]. However, based on the successes in antiviral drug discovery, 
recent years have seen an increasing interest in utilising antivirals not only for treatment, but 
also to prevent HIV transmission. Two general strategies are currently implemented for this 
purpose: (i) Treatment-as-prevention (TasP) intends to put individuals with an HIV diagnosis 
immediately on treatment, which essentially makes them non-contagious by decreasing the 
number of viruses they can expose to uninfected individuals [8]. (ii) Pre-exposure prophylaxis 
(PrEP) on the other hand prevents establishment of HIV infection after exposure [9, 10]. 

Currently, two oral PrEP options, the patent-expired two-compound combination Tru- 
vada, as well as the patent-protected two-compound combination Descovy are available. 
However, many more drugs are investigated for re-purposing [11, 12], or under de novo devel- 
opment [13], including topically-applied drugs, long-acting injectibles, as well as drug eluting 
implants [14, 15]. 

However, demonstrating clinical efficacy of novel PrEP compounds constitutes a formida- 
ble task. Clinical efficacy of PrEP is understood as the reduction in the number of infected 
individuals in a treatment- vs. a control arm of a clinical trial [9]. A major statistical problem 
arises from the fact that HIV transmission probabilities are extremely low (e.g. < 3% during 
unprotected sex [16]; far less for condom- or PrEP usage, and when potential transmitters take 
antiviral therapy). Hence, the number of evaluable data points (= infected individuals in a trial 
treatment and control arm) are extremely low and prone to chance events. Since the approval 
of Truvada-based PrEP, novel PrEP interventions have to be compared with Truvada, worsen- 
ing the statistical problem considerably [17]. E.g. the recent DISCOVER trial evaluating the 
efficacy of emtricitabine plus tenofovir alafenamide (Descovy) against Truvada was conducted 
over 8756 person-years [10], yielding as little as 22 evaluable data points (infections). 
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Statistically empowering such a study quickly exceeds organizational and monetary capacities. 
The statistical limitation has two consequences: (a) The determination of concentration-pro- 
phylaxis relations, threshold concentrations and contributions of transmitted and acquired 
drug-resistance cannot be rigorously deduced from clinical data, non-withstanding ethical 
concerns. (b) The hurdles to introduce next-generation PrEP regimen are immense: Trials 
consume huge monetary resources and require several years. This compromises the advance- 
ment of next-generation PrEP and likely affects its costs. It is therefore absolutely crucial to 
discern promising from less promising interventions a priori. 

Auxiliary tools based on integrative mathematical modelling may help to better understand 
the parameters contributing to clinical PrEP efficacy [18]. In particular, how drug dosing may 
alter the risk of acquiring HIV infection, depending on its timing and the magnitude of viral 
exposure. 

A key feature of HIV biology is that transmission is highly inefficient. For example, < 3% of 
unprotected sex acts between sero-discordant partners result in HIV infection [16]. Moreover, 
the number of genetically distinct founder viruses is extremely low [19]. This argues that sto- 
chastic processes play an important role during early infection, and that, in the majority of 
exposures, the virus becomes eliminated before it irreversibly infects the new host. Therefore, 
stochastic modelling and simulation approaches are used to estimate the efficacy of PrEP by 
integrating various host-, viral- and drug-specific parameters. For fixed drug concentrations, 
Monte-Carlo schemes, analytical, as well as probability generating ODE systems have been 
developed [20-22]. These approaches have been extended to include time-varying drug con- 
centrations by integrating pharmacokinetic characteristics, as well as realistic dosing schemes, 
but were restricted to particular drugs, drug classes or prophylactic schemes [23, 24]. Recently, 
a numerically exact Monte-Carlo approach was introduced that can be universally applied to 
study the effects of dosing, pharmacokinetics, drug adherence, timing and extend of viral 
exposure on the risk of HIV infection [11, 12]. Despite its advantage, the introduced stochastic 
simulation approach is still computationally prohibitive, in a sense that it would not allow to 
compute a PrEP efficacy from a history of drug dosing ‘on the fly’, e.g. to be useful in a health 
app or computer program that empowers PrEP users, akin to [25]. 

In this work, we derive a numerical method from an established viral dynamics model of 
HIV infection [26, 27] that overcomes aforementioned limitations. The method estimates the 
probability of viral extinction using a set of low-dimensional deterministic Ansatz functions 
that are solved with standard numerical solvers. The method allows to quantify PrEP efficacy 
within fractions of a second on a standard computer and is numerically exact up to the toler- 
ance level of an ODE-solver. This method fully integrates drug pharmacokinetics, which allows 
to estimate the influence of drug dosing, drug adherence, timing and magnitude of viral expo- 
sure on PrEP efficacy. 

We illustrate the method with the second-generation integrase inhibitor dolutegravir 
(DTG). Taking advantage of the outstanding performance of the developed method, we pre- 
sented several show cases to display their possible applications: By estimating (i) both pre- and 
post-exposure prophylaxis with different dosing of the drug, as well as timing- and extent of 
virus exposure, (ii) prophylactic protection profiles over a 6 month dosing history, as well as 
(iii) timing- and probability of viral establishment with different exposures. 


Methods 


The initial replication of HIV after exposure is highly stochastic. Typically, a low number of 
founder viruses is responsible for establishing infection and the transmission probability per 
sexual exposure is very low. In biology two types of stochastic noise are considered, roughly 
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categorized as internal- and external noise. Herein we focused on the internal noise, which 
assumes that the stochastic outcome of viral exposure can be explained by the order in which 
the considered reactions occur. For example, when a virus comes in close proximity of target 
cells, it can be either eliminated or infect a target cell and replicate. Prophylactic drugs shift the 
balance of these two events towards viral clearance. Stochastic dynamics of this type can be 
defined by a multivariate Poisson process, and the evolution of the state probabilities can be 
then described by the chemical master equation (CME). Conditioned on an initial state xo, the 
CME can be written down for any particular state x; as: 

d K 
nace =X; | X= xy) _ a a, (x; _ v,) . P(X, =X; — Vy | X= x») 

~ a,(x;) : P(X, = Xx; | X= Xo) 

with t > 0 and v; denoting the stoichiometry of reactions yielding state x;. Essentially, the 
CME as stated above, locally describes the flux of probability into- and out of the state x; with 
propensity functions a;. X, € N* denotes the state of the system (the combination of numbers 
of viral compartments, as well as drug particles), where s denotes the overall number of vari- 
ables. Because each variable can take any value in the natural numbers (e.g. 0, 1, 2, . . ., 00), the 
CME denotes an infinite dimensional system of ordinary differential equations, that is intrac- 
table (curse of dimensionality). 

In the context of estimating PrEP efficacy, we are however not interested in the entire state 
space but only interested in estimating the probability of extinction of all viral compartments, 
as outlined below. Thus, the key idea of our proposed method is to derive a low-dimensional 
system of ordinary differential equations, that allows to estimate only those probabilities rele- 
vant to estimating PrEP efficacy. 


Prophylactic efficacy 


The prophylactic efficacy is defined as the reduction of infection risk for a prophylactic regi- 
men S, compared to the infection probability in the absence of prophylaxis: 


PY, 8) 


9(Y,,S) =1 ~ PAY,, DB) 


(1) 
where P,(Y,,S) and P,(Y,, @) denote the infection probabilities in the presence and absence 
of a prophylactic regimen S for a given virus state (e.g. exposure) Y; at time t. We consider a 
prophylactic regimen to be a continuous function of drug concentrations. The state of the viral 
compartments Y; is defined as Y; = [V, T}, T2]", where V, T, and T> are the numbers of viruses, 
early stage infected T cells T, and late infected T cells T, respectively. For the absence of pro- 
phylaxis, P,(Y,, @), analytical solutions have been presented in [28]. For a prophylactic regi- 
men S they need to be determined numerically. 

The infection probability is the complement of the extinction probability Pz. We thus have 


P,(Y,,S) 7 1—P,(Y,,S) (2) 


where 
t 


0 V 
P,(Y,) :-=P[Y,=|0]| | y,=|T7, 
0 r, 


In words, the probability that all viral compartments will eventually go extinct, starting 
from state Y; at time ¢. Under the reasonable assumption of statistical independence, akin to 
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[28], we can define the extinction probability as: 
P,(Y,,8) = (P.(Y, = V,S))" : Pa, = TS) . (PY, = T.S\y? (3) 


where V, T,,, T, represent the unit vectors: 


1 0 0 
V=]0|, =]1|, F=|0 (4) 
0 0 1 


Thus, if P,(Y, = V,S),P,(Y, = T,,S) and P,(Y, = T,,S) can be determined, the prophy- 
lactic efficacy can also be calculated using Eqs (1)-(3). 

Below, we will present a method to compute the extinction probability numerically. First 
we will introduce a within-host HIV dynamics model, as well as a pharmacokinetic-pharmaco- 
dynamic model of the second-generation integrase dolutegravir, which serve as a common 
basis to derive and demonstrate the presented numerical method to compute HIV prophylac- 
tic efficacy. We will then introduce the proposed method (the Probability Generating System; 
PGS). A formal derivation of the method can be found in S1 Text. A pseudo-code for the 
method is found in $2 Text. To highlight the generality of the method, we derive the equations 
for the PGS for a different viral dynamics model that involves latently infected cell dynamics 
in S3 Text. 


HIV viral dynamics model 


We adapted the viral dynamics model from [26, 27]. We use this model, because it allows to 
mechanistically integrate the mechanisms of action of all approved drugs (and drug classes) 
[26] and because it allows to integrate both drug-specific in vitro and in vivo parameters [24]. 
In its most basic form, the considered viral replication cycle consists of free infectious viruses 
V, early infected T-cells (T,), and productively infected T-cells (T,). The dynamics can be 
defined by six reactions R, to Rg with propensities a,-a¢: 


R, : Clearance of free virus, V—-¥* a,=(CL+CL,-T,)-V (5) 
R, : Clearance of T,-cell, T,—>* a, = (dpc + 6;,) > T; (6) 
R, : Clearance of T,-cell, T, — * a, = 6, °T, (7) 

R, : Infection of a suscept. cell, V—-T, a,=f-T,-V (8) 
R, : Integration of viral DNA, T,—-T, a,(9)=k-T, (9) 
R, : Production of new virus, T,—>V+T, a,=N;,- Th, (10) 


where we assume that the integrase inhibitor dolutegravir (DTG) inhibits proviral genome 
integration (reaction Rs), with details outlined below. Moreover, we assume that a T-cell con- 
tinuously produces viruses (with reaction rate Rg) until it is cleared (continuous virus produc- 
tion model). The basic model is depicted in Fig 1. Utilized model parameters and their 
interpretations are given in Table 1. 
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Fig 1. Schematic of the utilized viral dynamics model. V, T,, T2 denote virus, early infected T-cells and productively 
infected T-cells respectively. Each reaction is denoted by its reaction propensity a, — dg. Briefly, a free virus can either 
be cleared (with reaction propensity a,), or infect a susceptible T cell with rate a, to yield an early infected cell T). 
These cells denote a state where the virus has penetrated the host cell, but has not yet integrated its proviral DNA into 
the host cell’s genome, thus not yet producing viral offspring. Early infected cells T, can either be cleared with rate a, 
or the proviral DNA irreversibly becomes integrated into the host cells DNA with rate as to yield a productively 
infected T-cell T,. T, cells start producing infectious progeny virus with rate ag, or they may get cleared by the immune 


system with rate a3. 





https://doi.org/10.1371/journal.pcbi.1009295.g001 


Pharmacokinetics and pharmacodynamics of dolutegravir 


Dolutegravir (DTG) is a second-generation integrase inhibitor which may potentially be used 
as prophylaxis against HIV. Moreover, we study it because of its similarity to carbotegravir, 
which is in clinical development as a next-generation PrEP compound. To evaluate the PrEP 
utility of DTG, we utilize a previously developed pharmacokinetic-pharmacodynamic model 


of the drug [11]. 


Pharmacokinetics of dolutegravir. We utilize the non-linear mixed effects pharmacoki- 
netic model introduced in [11]. In brief, a two-compartment model with first order absorption 
describes the plasma concentrations time profiles of dolutegravir (DTG) after oral drug 
administration. Individual parameters for a population of HIV-negative individuals were sam- 
pled from the distributions defined in [11] (Table 2 therein). The structural pharmacokinetic 


Table 1. Parameters for viral dynamics model. 















































Parameter Description Value Reference 
CL clearance rate of free virus by the immune system 2.3 29, 30] 
CLy clearance rate of the free virus during unsuccessful infection is ( a _ 1) Bp 26] 

Tt, level of uninfected T-cells in the disease-free state Ty = Ar/ by 

Prev probability of successful reverse transcription 0.5 31, 32] 
B lumped rate of infection of T-cells 8-107 33] 

Ar birth rate of uninfected T-cells 2-10° 34] 

or death rate of uninfected T-cells 0.02 35] 
Optic rate of intracellular destruction of pre-integration complex (PIC); 0.35 32, 36] 
Sr, rate of clearance of T,-cells 0.02 35] 

5x, rate of clearance of T,-cells 1 37] 

k rate by which T,-cells are transformed into T-cells 0.35 32] 

Nr rate of production of infectious progeny virus 670 26, 35] 








All parameters are in units [1/day], except for A [cells/day] and £ [1/(day - virus)] 





https://doi.org/10.1371/journal.pcbi.1009295.t001 
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model is given by the following set of ordinary differential equations (ODEs): 











d 
re — —k, ss Zi (11) 
d d k CL/F,; F,, F,, 
D= Z = a . z / bio , Z Q/ bio | vA +4 Q/ bio, Z, (12) 
dt dt Vd Pi Vid Foie Ve Fin Vol Pe 
d F,, F,, 
Ze = Q/ bio . z, Q/ bio, Zz (13) 
dt V./ Frio Vd Pes 


where Z, represents the amount of drug in the dosing compartment, and Z; denotes the DTG 
concentration in the peripheral compartment. D = Z, is the DTG concentration in the blood 
plasma, i.e. the value of interest. In a therapy, the value Z, increases whenever a dosing event 
T, occurs: Z, = Z,; + dose. 

Pharmacodynamics of DTG. Since DTG is an integrase inhibitor, it acts intracellularly 
by preventing the integration of proviral DNA. This effect can be translated into a decrease in 
propensity function a; by a factor 1 — np 


a;(t) = (1—mp(t))-k-T, (14) 


where 77p(t) denotes the direct effect of DTG at time ¢, which is modelled using the Emax- 
equation [38]: 
DI 


L) = 15 
m= ee (15) 


where D, is the drug concentration in the blood plasma at time t. ICs) represents the plasma 
drug concentration by which the activity of proviral integration is inhibited by 50%, and m 
denotes a hill coefficient. In this work ICs9 = 89 [nM] and m = 1.3 is used, which are values 
after protein adjustment [28] (free drug hypothesis). 


Low-dimensional deterministic Ansatz function 


In this section, we present a method to compute the extinction probability for the unit vectors 
P.(Y, = V, 8), P,(Y, = T,,S) and P,(Y, = T,,S), which enable the integration of arbitrary 
pharmacokinetic profiles resulting from some prophylactic regimen S. As noted before, this 
would enable calculating prophylactic efficacy for arbitrary drug/dosing regimen. 

Distribution of state transition events. The viral dynamic model illustrated in Fig 1 is 
interpreted as a continuous time, discrete state Markov process [39]. Therefore, the time when 
a particular state transition happens is exponentially distributed according to the reaction pro- 
pensities. In the viral replication cycle, if the initial state is Yo = [1, 0, 0] T ‘there are two possible 
next states: (i) the virus is cleared Y, = [0, 0, 0] a or a T-cell emerges Y; = [0, 1, ol’, where T is 
a random, exponentially distributed waiting time, until a reaction fires. The probability density 
function (PDF) for state transition V — T, can be derived as: 


Fyn, (*) = (1 — Fi, ()) - Sa, (2) 
= (1—(1—e™“*)) - ae") (16) 


— (ay +44 )x 
a,e 


where F, (x) is the cumulative probability of state transition V — © between time point 0 and 
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x, and f, (x) is the probability density function for transition V — T,. In words: The probabil- 
ity that V — T, occurs and V — @ has not occurred yet. Corresponding derivations hold for 
the process T, — T2 + V: 


fom) = (1-F,())-f,@) 
= (1- (1 e®*)) « (agen) (17) 


=a e (4346 )x 
6 


For the process T, — T>, the probability distribution f,, (x) is different since the values of as 
are time-dependent when an integrase inhibitor is applied (as in our example) that affects the 
reaction according to its pharmacokinetic-pharmacodynamic (PK-PD) properties and its dos- 
ing history. Using a Taylor approximation, the probability distribution for f,,_,,, (x) can be 
derived as (S1 Text): 


hao = = F,@)) 7,6) = aoe Ol a0 dt) (18) 


Probability Generating System (PGS). The precondition for the PGS method is that the 
concentration profile of the prophylactic regimen S must be known in advance. This can be 
achieved by solving the deterministic pharmacokinetic Eqs (11)-(13) for a particular dosing 
schedule using standard ODE-solvers. Given these pharmacokinetic profiles, the PGS method 
delivers an extinction probability-time profile. Mathematically, the method calculates 
P,(Y,,8). 

First, we derive the discrete form of this method. The discrete form is based on a time-dis- 
cretization of the underlying time-continuous Markov process into time steps At. From the 
constructed discrete-time Markov chain, we can compute the extinction probabilities for e.g. 


state V at time point t as follows: 


Py, = V) — PUY xj = 0 | ¥ = V+ 
PUY at = z, | i= V) “PEN gts = T,)+ (19) 
PUY gy = V | Y, = V) PT ai > V) 


=0 | Y, = V) denotes the probability that the virus 

= T, | Y, = V) is the probability that the state 
transition V — T, occurs in [f, t + At) and P,(Y,,,, = T;) is the probability that a T-cell that 
exists at time ¢ + At will eventually be eliminated. Finally, the probability that no state transi- 
tion occurs in [t, t + At) is given by 


with the following interpretations: P(Y,, ,, 


is eliminated in time span [t, t + At); P(Y,,., 


P(Y, 


t+At 


=V|¥,=V)=1=P0,4=0|¥,=V) PY, 


t+At 


=T,|Y,=V). 


The terms P(Y,,,, = 0 | Y, = V) and P(Y,,,, = T, | Y, = V) can be derived based on the 


calculations for the distribution of state transition events, as outlined in S1 Text. Using these 
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derivations, the the final expression for P,(Y, = V) is 

















“ a 
PAY,=V)= 1 1 —(a,+a4)At 
e( t ) a, +4, e ) 
a —(a,+a,)A a 
ah a 77 i, (1 —e (a, +a4) ‘) F aaa = T,) (20) 
ee MEE, PalYiias = V) 
Similarly, for P,(Y, = T,) we get: 
B a. 
P,(Y,=T,) = 2 (4 — prlmtas()ae 
(Y, = T,) nae ) 
a;(t) eas : 21) 
5 1— e (ay+a5(t))At) , PAY a7 ( 
a, ae) ) 2 ( tat ») 
A eran tag (At , PY ea, -_ T,) 
and for P,(Y, = T,): 
P,(¥, = T,) = —“3— (1 — enero) 
E\f¢ 2 7 
A —(agtae )A A 
ata (hee) PeFaw = 8) (22) 
P,(Yiea = V) + ectotot. BLY, = Ty) 


Having Eqs (20)-(22), we will take At — 0 to derive a set of continuous Ansatz functions. 
With this idea in mind, we differentiate Eqs (20)-(22) to derive the following set of ordinary 
differential equations (S1 Text): 


aP,(Y, = Vv) _ 4 
dt =a,-(P,(Y,=V)—1) 
+ ay (PAY, = V) — P,(Y, = T,)) 
dP,(Y, =T,) ; ie 
dt 2 (P,(Y, = T,) — 1) - 
+a,(t) - (P;(Y, = T,) — P,(Y, = T,)) 
rin ges Le _ 
dt a,-(P.(Y, = T,) — 1) +4, - (P;(Y, = T,) 


Given a;(t) and initial values, equations above can be solved by any ODE solver, as outlined 
in S2 Text. 
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Implementation and availability 


Pseudocodes for the method are described in detail in S2 Text. The method was implemented 
in Python 3.8, using SciPy 1.5.0, Numpy 1.18.5, Pandas 1.0.5 and matplotlib 3.2.2. Codes are 
available from https://github.com/KleistLab/PrEP. git. 


Algorithmic specifications 


Probability Generating System (PGS). The simulation time horizon must exceed the 
time horizon of interest. In this work, we extended the end time point for an extra t = 100 
hours so that the values in the target time interval are accurate. This value of extra time can be 
determined roughly from the half life of considered drugs, e.g. t = 7 - ty2, where ty/2 = 14.5 
hours denotes the half life of DTG. The above stated criteria guarantees that the drug concen- 
trations at T, + tT are < 1% of the trough concentrations. The method is solved backwards 
from an end time T, + tT to some start time T,. At T, + 7, we initialize the ODE with the values 
in the absence of drugs, e.g. P,(Y;, = V,S) =P,(V, 2), P,(Y7, = T,,S) =P,(T,, 2) and 
PCY, = T,,S) = P,(T,, @), which can be determined analytically [21]. We use sol- 
ve_ivp in SciPy [40] with the LSODA solver (linear multistep method) and default settings 
to solve the system of ODEs backwards in time. 

EXTRANDE. We implemented the exact stochastic simulation method EXTRANDE with 
configurations identical to [11]. 

Pharmacokinetics. Pharmacokinetic profiles for DTG were pre-computed from Eqs 
(11)-(13) using solve _ivp in SciPy [40] with default solver and default settings, for the 
respective prophylactic regimens S and D, was linearly interpolated. Using the pharmacoki- 
netic profiles, the time-dependent value of 7p(t), Eq (15) was determined and used in the 
respective algorithms. Depending on the analysis, we either simulated DTG pharmacokinetics 
for a representative individual, or by drawing pharmacokinetic parameters for 1000 virtual 
individuals from the parameter distributions defined in [11]. 


Simulation of pre- and post-exposure prophylaxis 


Single viral challenge. For single viral challenges the profiles of P,(Y, = V,S), P;(Y, = 
T,,S) and P,(Y, = T,, S) were computed using PGS (Eq (23)) with configurations outlined 
above. To compute extinction probabilities with arbitrary viral challenges we used Eq (3) and 
computed prophylactic efficacy from there using Eqs (1) and (2). 

Multiple viral challenges. Using the new methods, we can also compute the prophylactic 
efficacy following multiple viral challenges: Under the assumption of statistical independence 
[28], the extinction probability following multiple viral challenges is the product of extinction 
probabilities for single viral challenges at their corresponding time points. Le.:, 


9(Y, S) =—|-— PV 8) (24) 
m PA(Yay, @) 
with {t;}, i= 1, ..., 1 denotes a set of n viral exposures and 


PY 4.3) =1- []?™,) 


due to statistical independence, where P,(Y,,) is computed as described above. 
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If n viral exposures occur at the same time ft; = t, we have 
P,(Yy.5) = 1 — (P,(Y,))". 


To assess the validity of this assumption we also simulated multiple viral challenges with 
EXTRANDE and compared the results. 


Computation of density function of the extinction event 


So far, we computed the probability that viral extinction is eventually happening after viral 
exposure at some time f; P,(Y, = V,S). Le., we only care if extinction eventually occurs, 
regardless of when it happens. 

The introduced method can however also be used to estimate the probability that the 
extinction occurs in a specific time range. To solve for this probability, we can e.g. alter the ini- 
tial conditions of the PGS. In essence, we are interested in the extinction probability within a 
time range [fo, t.]. Since Eq (23) are solved backwards in time (S2 Text), we initialize the set of 
ODEs with the probability that the extinction occurs at the final time point f,, for example 
PY, =0/Y, = V) for the first ODE (Eq (23)). This probability is obviously 0. Consequently, 
the set of ODEs Eq (23) is initialized with values [0, 0, 0] and then solved backwards in time 
until fo. The derived probabilities have the following interpretation: they denote the probability 
that a single virus, single T, and single T>, transmitted at time point fp are cleared in the time 
range [fo, t.]. For example, if we want to determine the probability that virus is cleared t, = 10 
days after exposure to one infectious virus (at day fo) for some prophylactic regimen, we set 
the time span of the ODE set Eq (23) to [0, 10], initialize the ODE with [0, 0, 0]? at t, = 10 and 
solve the PGS backwards to fo = 0 to derive P(Y, =0,S | Y,, = V) in the first equation of the 
PGS, Eq (23). For a given pharmacokinetic (PK) profile D, this process has to be repeated for 
different values of t, to reconstruct the entire cumulative probability density function (CDF) 
with the correct PK profile. The probability density function (PDF) can be derived straight- 
forward from the CDF. Similarly, the procedure also allows to compute the corresponding 
probabilities for arbitrary numbers of initial viruses (statistical independence assumption), 

Eq (3). 


Results 
Model predictions in the absence of drugs 


To assess the validity of the proposed method and the utilized model, we initially performed 
some simulations without any drugs. 

Using the PGS, we calculated the probability of viral extinction when one virus 
P,(Y, = V, ), one early infected cell P,(Y, = T,, @) cell and one late infected cell 
P,(Y, = T,, @) were inoculated respectively. 

The respective extinction probabilities were 90.17%, 52.07% and 1.50%, respectively, which 
is identical to the analytical solutions (constant drug concentrations including drug absence) 
presented in [21]. 

We then computed the density function of the extinction event, as shown in Fig 2A using 
different inoculum sizes. Essentially, when a single virus reaches a replication-competent com- 
partment (red line in Fig 2A), if extinction occurs, it occurs very early after exposure, with the 
probability declining exponentially. When 10 viruses are inoculated (green line in Fig 2A), 
viral extinction is most likely to occur 14 hours after exposure. Overall extinction is less likely 
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Fig 2. Model simulations in the absence of drugs. A. Probability density function of the extinction event in absence of 
drug, for one initial virus and 10 initial viruses respectively. B. Relation between the expected inoculum size and the 
corresponding average infection probability per exposure, in the absence of drugs. The distribution of inoculum sizes is 
described in [24]. Average infection probabilities for homosexual vs. heterosexual transmission [16], are marked by dots. 








https://doi.org/10.1371/journal.pcbi.1009295.g002 


to occur with 10 viruses than with an inoculum of 1 virus. With both inocula, if viral extinction 
happens, it has to happen before day 5 post-exposure. 

Finally, we wanted to estimate the average infection probabilities Ey [P,(Y,,@)] = 
1—E, [P,(Y,, )] for different inoculum sizes and compare with published data, Fig 2B. To 
that end, we used viral load distributions P(VL = k) in potential transmitters, combined with 
the virus exposure model from [24] to computed average infection probabilities in case of het- 
erosexual vs. homosexual intercourse. 






































n 
2,,[P(¥,,2)] =>, | ¥,= | 0|,2 / P(VL =k) - P(V, = nlVL =k) 
n=0 k=0 
0 ————xXsX——’ 


probability to inoculate n viruses 


The simulations with the PGS in the absence of drugs (Fig 2B) indicate a clear increase of 
the infection risks with increasing virus exposure. Moreover, the predictions capture reported 
average per-exposure infection risks of 3% (homosexual intercourse) vs. 0.3% (heterosexual 
intercourse) [16] quite accurately. 

Next, we evaluated predictions, where we utilized actual pharmacokinetics of DTG to com- 
pute the prophylactic efficacy with the PGS method for different dosing regimen. 


Efficacy of PrEP-on-demand with DTG 


Using the proposed method, we computed the time course of the extinction probability 
P,(Y,,5) and the corresponding prophylactic efficacy y(Y,,S) for a 3-days once daily short- 
course oral 50mg DTG prophylaxis, that was either initiated shortly after viral exposure (post- 
exposure prophylaxis, PEP) or before virus exposure (pre-exposure prophylaxis, PrEP). Fig 3A 
shows the profiles of the extinction probability and Fig 3B depicts the corresponding prophy- 
lactic efficacy. In Fig 3, for illustration, we depict these quantities for Y, = V,Y,= T and 
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Fig 3. Prophylactic efficacy for 3-days once daily short-course oral 50mg DTG. The extinction probabilities were 
computed by PGS, for a representative individual with pharmacokinetic parameters k, = 2.24h'", VI Foio = 0.73 L, Q/Foio = 
0.0082 L/h, CL/Fpio = 0.85 L/h and V,./Fpio = 17.7 L. Observation began from two days before the first dose (the drug-doses 
are marked by arrows), until the 10th day after the first dose of DTG. The X-axis denotes the timing of viral exposure relative 
to the first dose, ie. negative values represent a viral exposure before the first dose of DTG (post-exposure prophylaxis, PEP), 
whereas positive values represent pre-exposure prophylaxis (PrEP) scenarios. A. DTG plasma concentration (dashed brown 
line) and the extinction probability profiles, with regards to one virus P,(V, S), one T;-cell P,(T, ,S) and one T-cell 
P,,(T,,S) are shown by solid red, green and blue lines. B. DTG plasma concentration (dashed brown line) and the 
corresponding prophylactic efficacies for one virus g(V, S), one T-cell g(T, ,S) and one T-cell g(T, , S) are shown by 
solid red, green and blue lines. 


https://doi.org/10.1371/journal.pcbi.1009295.g003 





Y, = T, individually. From these quantities, the extinction probability for an arbitrary initial 
state can be calculated based on Eq (3). 

From a computational point of view, PEP and PrEP are computed within the same execu- 
tion of the proposed methods: A pharmacokinetic trajectory (brown dashed line in Fig 3A and 
3B) is placed on the time axis and the methods are backwards propagated to some time before 
the first dose of the drug was given. Any time points before the first dose denote PEP, whereas 
PEP refers to time points after the first dose of the drug. For example, the value of the red line 
in Fig 3B at t = —2 (days) indicates that the efficacy of a 3 days 50mg oral DTG post-exposure 
prophylaxis, that was initiated 2 days after exposure to a single virus particle, is about 35%. If 
exposure occurred at t = 0 (coinciding with the first drug intake), the efficacy would be © 90%. 
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As described above, the new methods therefore compute a prophylactic efficacy profile for a 
time span of interest, i.e. every point on this time-efficacy curve represents the prophylactic 
efficacy p(Y,,S,,) conditioned that the viral exposure occurred at the indicated time point and 
the prophylactic regimen S was started at fy = 0. 

From a biological point of view, Fig 3B nicely highlights the role of the molecular target of 
the integrase inhibitor DTG within the viral replication cycle [28]: The drug is able to potently 
prevent infection if it emanates from a virus V or a T, cell (compartments preceding its molec- 
ular target process), but not from a T, (a compartment succeeding its molecular target 
process). 

Ona AMD R35 core with 3.6 Ghz and 16 GB RAM, PGS ran in a fraction of a second for the 
presented example in Fig 3. Notably, if stochastic simulation was performed (as in [11]), sev- 
eral thousand stochastic samples need to be generated to approximate the prophylactic efficacy 
¢(Y,,S) from the sample statistic for a single time point ¢. 


Comparison with stochastic simulation (EXTRANDE) 


Previously, an exact hybrid stochastic simulation method called EXTRANDE was introduced 
in [41]. We subsequently adapted the algorithm to estimate PrEP efficacy against HIV [11, 12]. 
Here, we used EXTRANDE to verify the accuracy of the proposed method. Fig 4 shows the 
predicted efficacy of a three days prophylaxis with either 2 or 50mg oral DTG started at t = 0 
using EXTRANDE vs. PGS. In contrast to Fig 3 (drug-centered evaluation), here we perform 


an exposure-centered evaluation to calculate 9(Y,, = V,S,). Le. the virus exposure occurs 
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Fig 4. Prophylactic efficacy for PrEP and PEP computed by EXTRANDE and the new methods. 2 mg (lines below) 
and 50 mg (lines above) DTG were ingested for three days, respectively. The extinction probabilities were computed by 
EXTRANDE and PGS, for a representative individual with pharmacokinetic parameters k, = 2.24h"’, Vp! Foio = 0.73 L, Q/ 
Fyig = 0.0082 L/h, CL/ Fig = 0.85 L/h and V,/ Fig = 17.7 L. The prophylactic efficacy was computed using Eq (1) for initial 
state Yo = [1, 0, 0]” (exposure to a single virus particle). The X-axis represents the timing of the first DTG dose relative to 
the virus challenge, which is marked by the arrow. EXTRANDE was run 10 000 times for each condition. The error bars 
denote the 95% confidence bounds for the ensemble estimate, computed using the Greenwood’s formula [42]. 


https://doi.org/10.1371/journal.pcbi.1009295.g004 
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to = 0 and the first dose is taken at time t; € {-23, -18, -12, —6, —3, -1, 2, 4, 6, 12, 18, 24} hours 
before/past the viral exposure. In words: the prophylactic efficacy if exposure to a single virus 
particle occurred at time ty = 0 and the 3-day prophylaxis was initiated at time t;. 

In Fig 4, we can see that the results of PGS form smooth lines, whereas the results of 
EXTRANDE fluctuate randomly around the results of PGS. From a biological standpoint we 
see that the prophylactic efficacy 50mg DTG is almost double compared to 2mg. Moreover, we 
see that the prophylactic efficacy deteriorates much faster for PEP than for PrEP. Le., if DTG is 
taken as PEP, it needs to be taken shortly after the exposure. For PrEP-on demand, the efficacy 
changes only marginally if PrEP is initiated within 24 hours prior to exposure. 

In terms of run time, EXTRANDE requires thousands of simulations to achieve statistically 
reliable meaningful results. In Fig 4, for each of the 12 time points we ran 10 000 simulations 
with EXTRANDE, which took about two hours for all points using multi-treading (12 treads). 
By contrast, using the proposed method, the values of all time points can be extracted from a 
single run, ie. p(Y,, = V, S,)=9(Y_, = V, S,,), with run times less than one second. 

Obviously, the pharmacokinetic profiles can be arbitrarily altered, which allows to assess 
the prophylactic efficacy of any regimen S, e.g. with regards to the drugs taken, their dose, the 
administration frequency and the timing of drug intake, as outlined above. 


PrEP efficacy for multiple viral challenges and different inoculum sizes 


Another interesting application of the proposed method is to assess the impact of the exposure 
on the prophylactic efficacy, e.g. to assess the sensitivity of p(Y,,S,,) with regards to Y;. We 
have shown previously by simulations [23] that the prophylactic efficacy depends on the inoc- 
ulum size (= how many viruses enter a replication-enabling environment). Also, in [24], we 
used stochastic simulation to assess prophylactic efficacy after multiple viral challenges. Here, 
we demonstrate how the proposed method can be used to address these questions. 

Multiple challenges. Multiple viral challenges can be computed straight-forward as exem- 
plified in the Methods section. In Table 2 we show the estimated prophylactic efficacy for dif- 
ferent PrEP regimen (either 3 or 7 doses of 2 or 50mg DTG started at ty = 0) with multiple 
challenges to a single virus particle, computed using both PGS and EXTRANDE. Foremost, as 


Table 2. Prophylactic efficacy in case of multiple viral challenges. 





























Dose # Doses Exposure time Prophylactic efficacy [%] Probability of Infection [%] 
EXTRANDE PGS P,(Yiq319) P,(Y;,),2) 

2mg 3 1, 24h 39.64+1.05 40.52 11.11 18.68 
3 1, 72h 27.7241.13 27.02 13.63 18.68 
3 1, 24, 72h 29.90+0.91 29.64 18.76 26.67 
7 1, 24, 72, 144h 44.31+0.73 44.92 18.66 33.87 

50mg 3 1, 24h 82.71+0.59 82.58 3.25 18.68 
3 1, 72h 72.21+0.74 72.27 5.18 18.68 
3 1, 24, 72h 73.33+0.60 73.79 6.99 26.67 
7 1, 24, 72, 144h 87.31+0.37 87.25 4.31 33.87 























2 or 50mg oral DTG was ingested for 3 and 7 days respectively starting at to = 0, and 2-4 viral exposures occurred at the time “Exposure time’ hours after DTG initiation. 


During each exposure one infectious virus entered a replication-enabling compartment. The corresponding prophylactic efficacy was computed by PGS and 
EXTRANDE, respectively. For each condition EXTRANDE was run 100 000 times. The 95% confidence bounds of the EXTRANDE estimate was computed using the 
Greenwoods formula [42]. The probability of infection after multiple viral challenges with- and without drug (P,(Y;,,,S) and P,(Y;,,,,@)) were computed with PGS and 
are depicted in the last 2 columns. Utilized pharmacokinetic parameters were k, = 2.24h'', V>/Foio = 0.73 L, Q/Fpio = 0.0082 L/h, CL/Fpio = 0.85 L/h and V,./Fpio = 17.7 L. 


https://doi.org/10.1371/journal.pcbi.1009295.t002 
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a sanity check, we can see that both methods yield congruent results for all tested conditions. 
We can also see that higher dose, as well as a longer time course (seven vs. three days) of DTG 
dosing improves the prophylactic efficacy, even if more viral challenges occur during the 
short-course prophylaxis. Also, we observe a interesting interplay between the number of 
exposures and their timing: For example, if two exposures occur at 1 and 24h after the first 
dose of DTG vs. three exposures at 1, 24 and 72h after DTG initiation, we see a decrease in effi- 
cacy. However, when we compare two exposures occur at 1 and 72h after the first dose of DTG 
vs. three exposures at 1, 24 and 72h after DTG initiation, we see a slight increase in efficacy. 
This has the following reason: In the presented example, the prophylactic efficacy (relative 
risk) after multiple challenges is given by: 


@(Yaa.S) =1- Pry S) (25) 
- P(¥,3,2) 
where {t;}, i= 1, ..., n denotes a set of n viral exposures. 


Now, if P,(Y;,;,@) increases faster with the number of exposures than P,(Y;,,,5), a sce- 
nario may arise, in which the prophylactic efficacy (= relative risk reduction) may be higher, 
although more exposures happened. Ie. the contextual information, ‘when did the exposure(s) 
occur relative to the drug dosing’ is relevant. For example, if many exposures happened at times 
of almost full protection, the exposed person would be better off than if only a few exposures 
happened at times of low protection. 

Inoculum size. Fig 5 shows how the profile of prophylactic efficacy is affected by the num- 
ber of inoculated viruses. This can be calculated from the solution of the PGS is a straight for- 
ward way, akin to Eq (3) (statistical independence assumption). In Fig 5, we observe that 
different inculum sizes lead to an (exponential) scaling of the prophylactic efficacy, and that 


PEP | | | PrEP Inoculum size 


Prophylactic efficacy [%] 





—2 0 2 4 6 8 10 
Time of viral exposure relative to the first dose [day] 


Fig 5. Profiles of prophylactic efficacy for different inoculum sizes. The experimental setup was chosen identical to 
Fig 3. To calculate the prophylactic efficacy profiles, an exponential scaling is applied to the solutions of the PGS (Eq 
(3)). The solutions are plugged into Eq (1). 


https://doi.org/10.1371/journal.pcbi.1009295.g005 
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the efficacy deteriorates, when large numbers of viruses are able to reach a replication-compe- 
tent compartment. 


Long-term prophylactic efficacy 


Because of its superior computational performance, the PGS can also be applied to estimate 
prophylactic efficacy over very long time scales for population pharmacokinetics (Pop-PK). I. 
e., typically pharmacokinetic variability is described by statistical models, such as non-linear 
mixed effects models (NLME), e.g. in [11]. If the pharmacokinetic characteristics of an individ- 
ual are not known, a Pop-PK model may still be used to accurately capture likely pharmacoki- 
netic profiles in an individual, given a dosing history. The PGS would then allow to predict the 
profile of prophylactic protection if the individual was exposed to virus at any time during the 
observation horizon. Note that this type of analysis is usually not feasible with stochastic meth- 
ods, due to computational demands (for each time point in the profile, several thousand sto- 
chastic simulations would be required). 

In Fig 6, we show the estimated long-term efficacy profile for a chronic, 6-month once- 
daily 50mg oral DTG regimen for different levels of adherence. For each adherence level the 
computation was conducted on 1000 virtual individuals sampled from the Pop-PK model to 
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Fig 6. Long-term prophylactic profile with different levels of adherence. N = 1000 virtual patients were sampled from the 
pharmacokinetic parameter distributions defined in Table 2 of [11]. 50 mg dose of DTG/ day was ingested in this six- 
month-long regimen with adherence level of 0.75, 0.5, 0.33 and 0.25. The red line depicts the median predicted prophylactic 
efficacy, whereas the dark- and light grey areas present the quartile range and the 2.5%-97.5% range respectively. The 
prophylactic efficacy was computed for Yo = [1, 0, 0]’. 


https://doi.org/10.1371/journal.pcbi.1009295.g006 
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capture inter-individual differences in drug pharmacokinetics. The red line and the grey 
ranges denote the median, interquartile and 2.5%-97.5% ranges of prophylactic efficacy under 
consideration of inter-individual differences in the pharmacokinetic profiles. With PGS it took 
about 24 min in total to compute the 6-month prophylactic profiles for 1000 virtual individuals 
and a given sequence of dosing events (determined by the adherence level), i.e. less than 1.5s 
for each individual on an AMD R5 core with 3.6Ghz and 16GB RAM (standard laptop). This 
computation could also be easily parallelized, which would reduce the run time considerably 
(the entire simulation took about 5 min on the same computer with 12 threads). 


Density function of the extinction event 


Using the approaches outlined in the Methods section, it is also possible to compute when the 
actual extinction event happens after exposure. This is highly useful in determining how long a 
prophylaxis on-demand should be given. 

Fig 7 shows the cumulative probability of extinction, as well as the density function of the 
extinction event, computed using the PGS for a 3days 50mg DTG regimen that was initiated at 
to = 0, coinciding with viral exposure. In these simulations, akin to the last example, we sample 
virtual individuals from the Pop-PK model. 

Fig 7A and 7B depict the cumulative-, as well as the density function of the extinction event 
after exposure to a single initial virus. Fig 7C and 7D show the corresponding distributions 
after exposure to 20 viruses. From the figures, we can see that the cumulative probability of 
viral extinction is much lower after exposure to 20- compared to one virus (Fig 7A vs. 7C). 
Moreover, we can see that, after exposure to a single virus, extinction most likely occurs shortly 
after exposure, e.g. within 1-2 days. In contrast, when 20 viruses are inoculated, extinction is 
most likely happening at day 4 after exposure, and that it is still likely that extinction may 
occur up to 10 days after exposure when a 3days 50mg DTG prophylaxis is applied. Moreover, 
extinction is less likely to happen: After exposure to 20 viruses, extinction occurs in 10 days 
with about 75% probability (median), compared to 98% after exposure to a single virus. 

In other words, our modelling highlights that the prophylactic efficacy depends on the mag- 
nitude of exposure. Moreover, the duration to eliminate the virus is prolonged when more 
virus becomes inoculated. Essentially, this suggests that large inocula, which may occur after 
blood transfusions, needle stick exposures, or tissue rupture during sexual contact may require 
longer duration of prophylaxis to prevent infection. 


Discussion 


The prophylactic efficacy of novel drug candidates against HIV is determined by the complex 
interplay of enzymatic, cellular, viral, immunological, pharmacological, as well as behavioural 
factors [18]. Some of these factors can be described in experimental surrogate systems [43]. 
However, their complex interplay, which determines clinical efficacy, can usually not be fully 
described using in vitro or ex vivo experiments. Moreover, animal models for HIV prophylaxis 
are suitable for proof-of-concept studies, but may still be confounded by inter-species differ- 
ences, as well as differences in the viruses used [44]. Lastly, while clinical studies of HIV pro- 
phylaxis are useful to assess the relevant efficacy endpoint, they do not allow to disentangle the 
complex interplay between the aforementioned factors, because of inherent limitations in the 
study design, sample sizes and the inability to measure the joint interplay of parameters that 
determine clinical efficacy. Therefore, it remains a formidable challenge to identify factors that 
alter prophylactic efficacy, or parameters that can be improved by the development of novel 
drugs or drug formulations for HIV prophylaxis. 
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Fig 7. Cumulative probability and probability density function of extinction event. N = 1000 virtual patients were 
sampled from the pharmacokinetic parameter distributions defined in Table 2 of [11]. 50 mg dose of DTG / day was 
ingested for three days, first dose was taken at fy = 0, coinciding with viral exposure. A. Cumulative extinction probability 
for one initial virus. B. Probability density function of extinction for one initial virus. C. Cumulative extinction probability 
for 20 initial viruses. D. Probability density function of extinction for 20 initial viruses. The red and green lines depict the 
respective median values, whereas the dark- and light grey areas present the quartile range and 2.5%-97.5% confidence 
range, taking inter-individual pharmacokinetic differences into account. 


https://doi.org/10.1371/journal.pcbi.1009295.g007 


Mathematical models [11, 12, 21-24, 29, 45-47] have become an essential tool to comple- 
ment our knowledge about prophylactic efficacy of antiviral compounds, as they are able to 
put the different parameters in context and test their relevance for determining prophylactic 
efficacy. Stochastic simulation methods are currently the gold standard for estimating prophy- 
lactic efficacy from these models [11, 12, 21-24, 29, 46, 47]. Essentially, to estimate prophylac- 
tic efficacy with stochastic simulation approaches, a large number of stochastic trajectories is 
sampled and subsequently classified into infection or extinction events to derive a sample sta- 
tistic of, for example, the probability of infection P,(Y,,S). More recently, the extra reaction 
algorithm for networks in dynamic environments (EXTRANDE) [41] was adapted in [11] to 
couple pharmacokinetics with intrinsically stochastic viral dynamics following exposure, and 
to accurately classify stochastic trajectories. Despite the advantages of stochastic approaches, 
the disadvantages are also very clear: To obtain meaningful statistics, many stochastic simula- 


tions need to be conducted to accurately determine the sample statistics P.(Y,,S). The latter 
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makes the stochastic methods expensive in terms of overall computational time, such that 
many important in silico experiments are infeasible. In particular, the scope of sensitivity anal- 
ysis with regards to aforementioned factors in integrated, multiscale-models is usually limited. 
In this work, we introduce a low-dimensional approach to estimate the prophylactic efficacy in 
considerably less time, in a single run. We envision that the approach can greatly expand the 
scope of analysis with regards to estimating prophylactic efficacy, by allowing to analyse the 
long-term effect of prophylaxis, as well as performing sensitivity analysis. 

Stochastic simulation methods sample trajectories of the whole system, where any state 
may arise during simulation. However, in the context of prophylaxis, one is only interested in 
the probabilities of extinction (and its complement, infection). Therefore, for each state, only 
those parts that contribute to the extinction event need to be considered. Because the probabil- 
ity of extinction for an arbitrary state can be expressed with the extinction probabilities of the 
respective unit vectors (Eq (3); statistical independence), only extinction probabilities for unit 
vectors need to be computed. This is the main idea behind the proposed low-dimensional 
Ansatz function. 

The discrete form of the proposed approach is based on the discrete-time Markov chain of 
the underlying virus dynamics model. The time was discretized into fixed time steps At for 
which the probability flux is computed (Eqs (20)-(22)). In this form, time-dependent func- 
tions, e.g. ds, are approximated by a zero-order Taylor approximation for each time step At. 
The continuous form, Eq (23) can easily be derived from the discrete form by replacing the 
constant time step At by an infinitesimally small step dt. The method is therefore based on a 
continuous-time Markov process. In our implementation of the method, we solved the phar- 
macokinetics beforehand and then wrapped the values of a; into a function that can be called 
directly from within the PGS’s ODEs, which are solved backwards in time. It is notable that a 
backward ODE solver must be used in the implementation, and the use and configuration of 
different ODE solvers will have an impact on the accuracy and efficacy of this method. 

The set of ordinary differential equations (Eq (23)) derived for the PGS is related to the Kol- 
mogorov backward equations [48], which was used in [22]. The major improvement of our 
work is to combine this Ansatz with the PK-PD profile of the prophylaxis (Eq (14)) so that the 
prophylactic efficacy can be computed with arbitrary prophylactic dosing schedules. Notably, 
instead of-, or in addition to PK-PD profiles, it is also possible to include other time-dependent 
functions that are relevant to the infection process. For example, any time-dependent, ie. 
adaptive immune response could in principle be considered. 

In Fig 4, we compared the proposed method with EXTRANDE, which delivers consistent, 
but less accurate results. Likewise, Table 2 shows the consistency of prediction with PGS and 
EXTRANDE, regarding the prophylactic efficacy for different DTG regimens with multiple 
viral exposures. 

Using the proposed method, we can also perform analyses that are computationally infeasi- 
ble with stochastic simulations. As shown in Fig 6, the long-term prophylactic efficacy profiles 
were computed using PGS with four different values of adherence and a virtual patient cohort 
of 1000 individuals. This type of analysis allows to quantify sensitivity with regards to adher- 
ence and inter-individual pharmacokinetic variability. The analysis showed that DTG can pro- 
tect highly adherent individuals from acquiring HIV infection and that inter-individual 
differences are most strongly affecting prophylaxis at times of inconsistent use of the prophy- 
laxis. Moreover, when several consecutive pills are missed, prophylactic efficacy may drop 
below 50%. Interestingly, the long-term prophylactic efficacy computation took less than 2 sec- 
onds for one virtual patient. To compute the corresponding profile with EXTRANDE, for 
example if an efficacy estimate is to be computed for every hour using 5000 simulations, a total 
of 24 x 365 x 5000 = 44 million simulations for a single adherence level would need to be 
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conducted. Besides computational time, power consumption (and possible carbon imprint) 
could therefore be considerably reduced using the proposed methods. 

Another possible application of the proposed methods is the possibility to estimate the time 
of the extinction event, Fig 7. As a showcase of a sensitivity analysis, we estimated, for a 3-days 
50mg DTG regimen, the probability density function of the extinction event, when a single, 
versus 20 viruses were initially reaching a replication-competent physiological compartment. 
The analysis showed that the time to viral extinction is increased for larger inoculum sizes 
(more viruses). This analysis thus highlights the complex interplay between viral exposure and 
prophylaxis that can be analysed with the proposed methods and be used to optimize HIV 
prophylaxis. With regards to data, interestingly, in vaginal SHIV viral challenge models, which 
are typically conducted with large inoculum sizes, late viral breakthrough has been observed 
[49, 50]. 

The presented method has been derived for the model depicted in Fig 1 and need to be 
adjusted if other viral dynamics models were used. In S3 Text, we derive the method for an 
extended viral dynamics model that additionally considers infected T cells, which may turn 
dormant, e.g. become latently infected. It is well established that these latent reservoirs are a 
major obstacle to the elimination of HIV during therapy [51]. While these reservoirs are estab- 
lished early in infection, it is unclear whether they alter prophylactic efficacy. To test this 
hypothesis, we used the proposed methods for the extended viral dynamic model (S3 Text). 
When comparing the results with those from the simpler model, we however found that the 
impact of the reservoirs on prophylactic efficacy was negligible. 

In summary, we propose a novel method that can estimate the efficacy of arbitrary prophy- 
lactic regimen and viral exposures within seconds. The method allows to integrate individual 
PK/PD profiles and viral dynamics into a single framework and it is more exact than state-of- 
art hybrid stochastic simulation schemes, like EXTRANDE. We envision that the new method 
can be applied in many circumstances, in which the stochastic simulation is computationally 
infeasible, such as parameter sensitivity analysis or long-term efficacy estimation. To this end, 
the proposed method may even be suitable as part of an App, which may help PrEP users to 
monitor and plan their PrEP regimen. Moreover, the general schemes may be adapted to 
study related biomedical questions, like prophylactic efficacy in other pathogens or vaccine 
efficacy. 


Supporting information 


$1 Text. The supplementary text contains a supplementary derivations for proposed 
method (distribution of state transition and PGS). 
(PDF) 


$2 Text. The supplementary text contains a complete pseudo-code and the implementation 
details for the Probability Generating System (PGS). 
(PDF) 


$3 Text. The supplementary text entails the derivation of the equations of PGS for a 
extended viral dynamic model that contains latently infected viral reservoirs. 
(PDF) 
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